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Abstract 

We investigate a Maxwell model of inelastic granular mixture under the influence of a stochastic 
driving and obtain its steady state properties in the context of classical kinetic theory. The model 
is studied analytically by computing the moments up to the eighth order and approximating the 
distributions by means of a Sonine polynomial expansion method. The main findings concern the 
existence of two different granular temperatures, one for each species, and the characterization 
of the distribution functions, whose tails are in general more populated than those of an elastic 
system. These analytical results are tested against Monte Carlo numerical simulations of the model 
and are in general in good agreement. The simulations, however, reveal the presence of pronounced 
non-gaussian tails in the case of an infinite temperature bath, which are not well reproduced by 
the Sonine method. 

PACS numbers: 02.50.Ey, 05.20.Dd, 81.05.Rm 
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I. INTRODUCTION 



Granular materials, a term coined to classify assemblies of macroscopic dissipative objects, 
are ubiquitous in nature and play a major role in many industrial and technological processes. 

Interestingly, when a rarefied granular system is vibrated some of its properties are similar 
to those of molecular fluids, while others are unique and have no counterpart in ordinary 
fluids ||T| . A spectacular manifestation of this difference can be observed in a driven mixture 
of granular particles: if one measures the average kinetic energy per particle, proportional to 
the so called granular temperature, one finds the surprising result that each species reaches a 
different value. Such a feature, observed in a recent experiment 0], is in sharp contrast with 
the experience with other states of matter. At a more fundamental level such a behavior 
is in conflict with the Zeroth Law of thermodynamics. This principle states that when two 
systems globally isolated are brought into thermal contact they exchange energy until they 
reach a stationary state of mutual equilibrium, characterized by the same value of their 
temperatures. A corollary to such a principle is the statement that the thermal equilibrium 
between systems A and B, i.e. = T B , and between A and C (T4 = T c ) implies the 
thermal equilibrium between B and C ||. 

On the contrary, when two granular system A and B subject to an external driving 
force exchange energy they may reach in general a mutual equilibrium, characterized by two 
constant but different granular temperatures and Tg. In addition even when A and C 
are in equilibrium at the same temperature Tc and B and C also have the temperature Tc 
one cannot conclude that A and B would be in mutual equilibrium at the same temperature. 
In other words, one of the most useful properties of the temperature, i.e. the independence 
from the thermal substance is lost when one deals with granular materials. 

In the present paper we shall investigate the properties of a simple model of two com- 
ponent granular mixture. The motivation of our study relies in the fact that granular 
materials being mesoscopic objects are often constituted by assemblies of grains of different 
sizes and/or different physical and mechanical properties. The study of granular mixtures 
has attracted so far the attention both of theoreticians M and of experimentalists [0. In 
particular Garzo and Dufty have studied the evolution of a mixture of inelastic hard spheres 
in the absence of external driving forces, a process termed free cooling because associated 
with a decrease of the average kinetic energy of the system, i.e. of its granular temperature. 
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During the cooling of a mixture, which can be homogeneous or not, according to the presence 
of spatial density and velocity gradients, each species may have different granular temper- 
atures, although these may result asymptotically proportional, i.e. they might decrease at 
the same rate. 

Menon and Feitosa instead studied several mixtures of vibrated inelastic grains and re- 
ported their failure to reach the same granular temperature. 

In the present paper our interest will be concentrated on the statistically stationary state 
obtained by applying an energy feeding mechanism represented by a stochastic driving force. 

The most widely used model of granular materials is, perhaps, the inelastic hard sphere 
model (IHS) |J, which assumes the grains to be rigid and the collisions to be binary, instan- 
taneous and momentum conserving. The dissipative nature of the collisions is accounted 
for by values less than 1 of the so called restitution coefficient r. Even such an idealized 
model represents a hard problem to the theorist and one has to rely on numerical methods, 
namely Molecular Dynamics or Event Driven Simulation, or to resort to suitable truncation 
schemes of the hierarchy of equations for the distribution functions. One of these schemes 
is represented by the Boltzmann equation based on the molecular chaos hypothesis, which 
allows to study the evolution of the one-particle distribution function. Its generalization to 
the two component mixture of inelastic hard spheres has been recently considered by Garzo 
and Dufty. 

Our treatment will depart from previous studies, because we have chosen an even simpler 
approach based on the so called gases of inelastic (pseudo)-Maxwell molecules. These gases 
are natural extensions to the inelastic case of the models of Maxwell molecules ||, where the 
collision rate is independent of the relative velocity of the particles. Such a feature greatly 
simplifies both the analytical structure of the Boltzmann equation and the numerical imple- 
mentation of the algorithm simulating the gas dynamics. Although the constant collision 
rate is somehow unrealistic one may hope to be able to capture some salient features of gran- 
ular mixtures, and in particular to reach a better understanding of their global behavior, 
because the model lends itself to analytical studies. 

This type of approach to granular gases had recently a surge of activity since the work 
of Ben-Nairn and Krapivsky 0, |8j, our group |5], 110, O, Ernst and coworkers and 
Cercignani and collaborators |TB|, |TJJ, and is providing a series of new important results 
concerning the energy behavior and the anomalous velocity statistics of granular systems. 
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The paper is organized as follows: in Section 2 we define the model and deduce the 
associated Boltzmann equations governing the evolution of the velocity distribution functions 
of the two species and set up their moment expansion; in section 3 we shall the exact values of 
the stationary granular temperatures; in section 4 we consider the moments up to the eighth 
order and compute the distribution functions by means of the so called Sonine expansion. 
In section 5 we simulate on a computer the dynamics of the inelastic mixture and compute 
numerically the distribution functions and compare these with the theoretical predictions. 
Finally in 6 we present our conclusions. 



II. DEFINITION OF THE MODEL 



In the following we shall consider a mixture of Maxwell inelastic molecules subject to a 
random external driving. Let us consider an assembly of N± particles of species 1 and N 2 
particles of species 2. For the sake of simplicity we assume the velocities v f, with a = 1, 2, to 
be scalar quantities. The two species may have different masses, mi and m 2 and/or constant 
restitution coefficients which depend on the nature of the colliding grains but not on their 
velocities, i.e. rn,r 2 2 and ryi = Ti\. 

The mixture evolves according to the following set of stochastic equations: 



m i ^ = F i + f i + Z i (t) (1) 

where the total force acting on particle % is made of three contributions: the impulsive 
force, Fi, due to mutual collisions, the velocity dependent force, = —Tvi, due to the 
friction of the particles with the surroundings and the stochastic force, £j, due to an external 
random drive. Since we are interested in the rapid granular flow regime, we model the 
collisions as instantaneous binary events, similar to those occurring in a hard sphere system. 

The presence of the frictional, velocity dependent term in addition to the random forcing 



T5L not only mimics the presence of friction of the particles with the container, but also 



is motivated by the idea of preventing the energy of a driven elastic system (7 — > 1), to 
increase indefinitely. 

Let us observe that in the absence of collisions the velocity changes are described by the 
following Ornstein-Uhlenbeck process: 
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m a d t v i {t) = -Tv i {t)+^{t) (2) 

where the stochastic acceleration term is assumed to have a white spectrum with zero 
mean: 



m*)) = o (3) 

and variance: 



<e?(f)#(0> = 2i>^M(* - o ( 4 ) 

By redefining the bath constants T a = ^- and D Q = it is straightforward to obtain 
the probability distributions of the velocity of each species, P a (v,t). In fact, the Fokker- 
Planck equations associated with the process (0): 



d t P a (v, t) = T a (d v vP a (v, t)) + D a (d 2 v P a (v, t)) (5) 
possess the following stationary distribution functions: 



where T b = S represents the temperature of the heath-bath, that we fix to be the same 



for the two species [JT6 



In order to represent the effect of the collisions on the evolution of the system we assume 
that the velocities change instantaneously according to the rules: 



v't = v? - [1 + r aP ] ^ « - v?) (7a) 
m a + 171/3 J 

v'f = vl + [1 + r aP ]^^{v? - vf) (7b) 
J J m a + mp J 

where the primed quantities are the post-collisional velocities and the primed are the 
velocities before the collisions. A finite fraction of the kinetic energy of each pair is dissipated 



during a collision. Between collisions the velocities perform a random walk due to the action 
of the heat-bath. The typical time r c between particle-particle collision is assumed to be 
large compared to the time between random kicks. On the other hand the typical time 
scales associated with the bath are r&i = = jr- and r^ 2 = ^ = When r — > we must 
also take -D — > in performing the elastic limit, otherwise the kinetic energy would diverge 
asymptotically. In fact in section V we shall discuss the situation of inelastic particles with 
vanishing friction, already considered in refs. Q jl3fl . 



The evolution equations for the probability densities of finding particles of species a with 
velocity v at time t for the system subject both to external forcing and to collisions are 
simply obtained by adding the two effects: 



dtP a (y,t)=r a (d v vP a (y,t)) + D a (d 2 v P a (v,t)) + -Q a (Pi,P 2 ) (8) 

where the collision integrals Q a consist of a negative loss term and a positive gain term 
respectively: 

Q l (P h P 2 )=- Pl ( v ,t) + -^— [duP 1 (u,t)P 1 p^-(l-r n )^ \ 



+ v i pjmi + m 2 duPi ^ t) p 2 _rm \rm 



(9a) 



1 + r 12 m 2 J V 1 + r 12 

J- + r 22 J \ -L + ^22 / 

1 + r i2 mi J \ 1 + r i2 



(9b) 



where p = N\/(Nx + iV 2 ). In writing eqs. @ we have assumed that the collisions occur 
instantaneously and that collisions involving more than two particles simultaneously can be 
disregarded. Moreover, all pairs are allowed to exchange impulse regardless of their mutual 
separation. In this sense we are dealing with a mean field model. In order to proceed 
further it is convenient to take Fourier transforms of eqs. (|9|) and employ the method of 
characteristic functions [Ol defined as: 



P a (k,t) = I dve ihv P a (v,t) (10) 
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The resulting equations read: 



dA(k,t) = -D 1 k 2 P 1 {h,t) - r^ACM) - -[A(M) -pA(tiiM)A((i -m)Kt) 

-(l-p)A(TteM)A(l-7ia)M)] 

(11a) 

d t P 2 (k,t) = -D 2 k 2 P 2 (k,t)-T 2 kd k P 2 (k,t) - -[P 2 (k,t) - (l-p)P 2 (i 22 k,t)P 2 ((l- l22 )k,t) 

-pA(72iM)A((l-72i)M)l 

(lib) 



with 



7«tf = (12a) 
712 = 11-^(1-712)] (12b) 

72i = [1 - T^ T ( 1 " (12c) 

with ( = mi/m 2 

From the mathematical point of view the driven case is very different from the cooling 
case. As we have seen in the latter case the tails are originated by the presence of a singular 
point at the origin in the equation for the characteristic function [0. The singularity is 



due to the conspiracy of the constant cooling rate and of the scaling form of the evolution 
equation. Such a singularity is removed in the driven case, thus high velocity tails cease to 
exist and the distribution will be much more well behaved. 

The mathematical structure of eqs. ( |TT1) is particularly simple and in fact there exists a 
standard method of solution. It consists in expanding the Fourier transform P a (k,t) of the 
distributions P a (v,t) in a Taylor series around the origin k = 0: 

Pa(k,t) = J2 { ---r^(t) (is) 

n=0 

and substituting (|T3|) into eqs. (|TT|) . 

Equating like powers of k we obtain a hierarchy of equations for the /i"(t) which can 
be solved by a straightforward iterative method. At this stage one can appreciate the 



mathematical convenience of the Maxwell model. In fact, the coefficients of the Taylor 
series represent the moments of the velocity distributions 

/oo 
dvv n P a (v,t) (14) 
-oo 

Since the evaluation of the moments of a given order requires only the knowledge of the 
moments of lower order one can proceed without excessive difficulty to any desired order. In 
practice, we carried on our calculation up to the eighth moment, assuming that the initial 
distributions were even, so that the odd moments vanish. In order to render the reading of 
the paper more expeditious we shall report the equations determining the stationary value 
of the moments in the Appendix. 



III. TWO TEMPERATURE BEHAVIOR 

In order to determine the granular temperatures we equate the coefficients of order k 2 in 
eqs. ( |TTD and obtain the governing equations for the second moments: 

r c d tl if = {p[27u(7n " 1)] " (1 - P)(l " 7i 2 ) - 2r c r 1 }/4 1) 

(15a) 

+ (l-p)(l-7i 2 ) 2 /4 2) + 2r cJ D 1 
tAv? ={(1-P) [2722(722 - 1)] - p(l - 7 2 2 i) " 2r c r 2 }/4 2) 

(15b) 

+ p(l-72i) 2 /4 2) + 2r CJ D 2 
The r.h.s. of eqs. ([15]) represent the balance between the energy dissipation due to 
inelastic collisions and friction and the energy input due to the bath. We define the global and 
the partial granular temperatures respectively as T g = pTi + (l—p)T 2 and T a = \m a < v 2 >, 
where the average is performed over the noise £j. Since the energy dissipation and the energy 
supply mechanisms compete, the system under the influence of a stochastic white noise 
driving achieves asymptotically a statistical steady state. Notice that eqs. flilf) feature only 
the second moments of the velocity distributions, so that the solution is straightforward. 
Such a state of affairs should be contrasted with the analogue problem of determining the 
partial temperatures in Boltzmann models Q. 

Let us start analyzing the behavior of the Maxwell gas in the one component limit limit 
p — > 1. The granular temperature approaches its stationary value Ti exponentially: 
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Ttit) = T!(0)e-^ + 7i(oo)[l - e"*] 



(16) 



where the constant r represents a combination of the two characteristic times of the 
process given by: 

1 _ 7ll (l - 711) | r 
r r c mi 

which shows that TJ, is an upper bound to the granular temperature. We also obtain a 
simple relation between the temperature of the bath and the granular temperature T\. 



T ^ = m ^ — Sfa^i) ^ Y = Tb (18) 



On the other hand, when the two components are not identical eqs. fllq) show that the 
equilibrium macro-state is specified by two different partial granular temperatures, both 
proportional to the heath bath temperature. Hence, the temperature ratio is independent 
from the driving intensity D as one can see from the formula: 



Tl _ 1 (1 - p) [2 722 (1 - 722)] + p(l - 7 2 2 i) + 2£ + (1 - P)C 2 (1 - 712) 2 

(19) 



T 2 C p[2 7 n(l-7n)] + (l-p)(l-7i 2 2) + 2^+K- 2 (l-72i) 
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Formula in eq. (O) illustrates the two temperature behavior of an inelastic mixture 
subject to external driving. Notice that the temperature ratio in the driven case is different 
from the corresponding quantity in the cooling undriven case, for the same model system. 
In the undriven case we found that the homogeneous cooling state was characterized by two 
different exponentially decreasing temperatures, but whose ratio was constant. However, 
no simple relation exists between the ratio relative to the two cases, on account of the fact 
that the energy exchanges involved are rather different. Thus, in the presence of a heath 
bath the inelastic mixture displays the two temperature behavior already reported in the 



free cooling case 0, [18| and in experiments 0. This feature seems to be a general property 
of inelastic systems. In fig. (|I|) we display the temperature ratio as a function of the mass 
ratio ( for two different values of the inelasticity. Notice that the temperature of the heavier 
component is lower than the one of the lighter species. 
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FIG. 1: Temperature ratio as a function of the mass ratio for different choices of the inelasticity 
parameter, for p = 0.5 and rn = r22 = ^12 = r = 0.8 (solid line) and r = 0.5 (dashed-line) 

1.3 1 1 1 1 1 1 1 1 1 1 
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. .. r n =0.9 
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FIG. 2: Temperature ratio as a function of the asymmetry Tiijr\\ for p = 0.5, £ = 1 and different 
choices of the inelasticity parameter m, from top to bottom this is respectively 0.99, 0.9, 0.4 (line) 

In fig. the ratio of the temperatures of the two species is plotted in the case of an 
asymmetry in the restitution coefficients parametrized by the form: r2i = rn — x with 
r i2 = ( r 22 + r n)/2, for p = 1/2, identical masses and three different values of the coefficient 
r n as shown in figure. One sees that the variation of |f is much smaller than the corre- 
sponding variation with respect to the mass asymmetry shown in fig. [I] in agreement with 
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the experimental observation of reference . It seems reasonable to conclude that the mass 
asymmetry is the larger source of temperature difference between the two components. 



IV. VELOCITY DISTRIBUTION FUNCTIONS 

An interesting aspect of granular systems concerns the nature of the single particle ve- 
locity distributions. The inelasticity, in fact, causes marked departures of P a (v,t) from the 
Gaussian form which characterizes gases at thermal equilibrium. In undriven gases these 
deviations are particularly pronounced and one observes inverse power law high-velocity 



tails both in gases of pseudo-Maxwell molecules [Jg] and in IHS ||19|| . In the driven case, 
i.e. in systems subject to Gaussian white noise forcing (similar to that represented by eq. 
(0) with T = 0) exponential tails of the form exp(— v 3 ^ 2 ) have been predicted theoretically 
in inelastic hard-sphere models [|19j and tested by direct simulation Monte Carlo of the 
Enskog-Boltzmann equation |p20[. Have these non Gaussian tails a counterpart in Maxwell 



models? Ben-Nairn and Krapivsky [7| on the basis of a re-summation of the moment ex- 
pansion concluded that the scalar Maxwell models with vanishing viscosity (T = 0) should 
display Gaussian-like tails. However, this prediction is in contrast with the argument, em- 
ployed by Ernst and van Noije in the case of IHS, which consists in estimating the tails of 
the distribution by linearizing the master equation (|9|) by neglecting the gain term. This 
assumption simplifies the analysis and allows us to reach the conclusion that the velocity 
distribution for large v should vanish as: 



lim P(v) oc exp (—v/vo) (20) 

v—>oo 

with Vq = Dt c . Clearly such a result is in sharp contrast with the result of ref. and 
seems to indicate that the Sonine expansion does not reproduce faithfully the high- velocity 
tails in the case of Maxwell models with vanishing viscosity. The test of the limit ( p0[ ) will 
be shown in the section 6, where we illustrate the results of our numerical simulations. 

On the other hand, the same kind of asymptotic analysis sketched above, allows us to 
conclude that the presence of a viscous damping is the redeeming feature which renders 
convergent the Sonine expansion and the associated Gaussian tails. In fact, with a finite 
value of T the asymptotic solution is of the form: 
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lim P(v) oc exp (-Cv 2 ) (21) 

v — >oo 

We shall test such a prediction in the remaining part of this section and study the velocity 
distributions of the individual species when T 7^ by constructing the solution to the master 
equation using the Sonine polynomial expansion method, one of the traditional approaches 



to the solution of the Boltzmann equation |2l 



We shall also investigate whether the two partial distributions can be cast into the same 
functional form upon re-scaling the velocities with respect to the partial granular thermal 
velocity, in other words if it is possible to have a data collapse for the two distributions. 

We shall first obtain the steady state values of the first eight moments as illustrated in the 
Appendix and then compute the approximate form of the distribution functions by assuming 
that these are Gaussians multiplied by a linear combination of Sonine polynomials. 

Let us begin by writing the following Sonine expansion of the distribution functions: 

Uc) = ^e- c2 [l + J2<S n (c 2 )} (22) 
where f a (c) is the re-scaled distribution defined by: 



f a (c) = y%qP a (v). (23) 

and c 2 = v 2 /2fi2- The expansion gives the distributions in terms of the coefficients a" 
of the Sonine polynomials S n (c 2 ). In practice, one approximates the series (]22]) with a 
finite number of terms. Since the leading term is the Maxwellian, the closer the system 
to the elastic limit, the less term suffice to describe the state. The expression of the first 
polynomials is: 
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S (c 2 ) = 1 (24a) 
S 1 (c 2 ) = l-c 2 (24b) 
S 2 (c 2 ) = | - \c 2 + \c* (24c) 

In order to obtain the first m values a^, we need to compute the re-scaled moments (c n ) a 
of the distribution functions up to order 2m. These moments are evaluated in Appendix by 
means of a straightforward iterative method. At the end, knowing the re-scaled moments, 
one obtains the following relation for the coefficients: 



a n = (25) 



Eq. (|25|) can be proved by imposing the consistency condition: 



(c n ) a = / dcc n f a (c) (26) 

J — oo 

in conjunction with the orthogonality property of the Sonine polynomials: 

/OO 1 
—e- c2 S n (c 2 )S m (c 2 ) = AU m ,n (27) 
-oo V 71 " 

where M n is a normalization constant. Notice that in order to obtain our results we 
have not assumed weak inelasticity, therefore these hold for any value of the restitution 
coefficients. 

The coefficients up to the fourth order in terms of the re-scaled moments read: 
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< = 



^ = [l-4(c 2 ) a + -(c 4 ) ( 



[l-6(c 2 ) a + 4(c 4 ) Q --(c 6 ) 



32 
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[l-8(c 2 ) Q + 8(c 4 ) Q --(c b ) Q + — (c 8 ) a ] 
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FIG. 3: Second Coefficient of the Sonine expansion 02, for each component as a function of the 
mass ratio £, mi = 1 and for p = 0.5 and r\\ = T22 = r\i = 0.5. We have also fixed p = r&i = 200 
and r c = 25. 

In fig. (0-0) we display the behavior of the Sonine coefficients for both components in 
the case of equal restitution coefficients r = 0.5 as a function of the mass ratio. 

In figs. (0-§|) we illustrate the variation of the Sonine coefficients for the two components 
as a function of the inelasticity for two different values of the mass ratio: ( = 1 and ( = 2. 
Notice that the coefficients are monotonic functions of the inelasticity as already noticed in 
pure systems |T9|| . 

In fig - (0) we show the distribution functions for the heated system with non vanishing 
viscosity. The tails become fatter with increasing order of the approximation, i.e. the high 
energy tails are overpopulated. Moreover, one sees that when p — 1/2 and the restitution 
coefficients are all equal the species with the larger tails is the lighter. On the other hand, 
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FIG. 4: Third Coefficient of the Sonine expansion 03, for each component as a function of the mass 
ratio (. The remaining parameters as in fig. || 
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FIG. 5: Fourth Coefficient of the Sonine expansion 04, for each component as a function of the 
mass ratio £. The remaining parameters as in fig. ^ 

for a system with the same masses, but different restitution coefficients, the more elastic 
species displays the larger tails. We also show the numerical data obtained by simulating the 
the dynamics. The agreement is quite satisfactory and validates the approximation method 
employed. 

On the other hand, it is also evident, that the two distributions fail to collapse one over 
the other after the rescaling of the velocities. This fact is consistent with the different values 
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FIG. 6: Second Coefficient of the Sonine expansion 02, for each component as a function of the 
inelasticity r = r%i = r22 = T2, f° r P = 0.5 and mi = 1 and £ = 2 and £ = 1. The remaining 
parameters as in fig. ||. 
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FIG. 7: Third Coefficient of the Sonine expansion 03, for each component as a function of the 
inelasticity r = m = T22 = r i2, for p = 0.5 and mi = 1 and £ = 2 and £ = l.The remaining 
parameters as in fig. ^. 

assumed by the coefficients a* and a^. However, this effect is rather small and can be 
appreciated only by studying the high velocity region of the distribution functions. 
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FIG. 8: Fourth Coefficient of the Sonine expansion a4, for each component as a function of the 
inelasticity r = r\\ = ^22 = r i2 ; for p = 0.5 and m\ = 1 and £ = 2 and £ = 1. The remaining 
parameters as in fig. 0. 
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FIG. 9: Rescaled distribution functions for the two components for a model system characterized 
by mi = 1, 771.2 = 0.5 and equal inelasticities r = 0.5.2], = 1 and the remaining parameters as in 

fig.| 

V. NUMERICAL SIMULATIONS 



To investigate the validity of the previous results and in particular to test the convergence 
of the Sonine expansion in different situations we shall present in this section numerical 
results obtained by simulating an ensemble of iV particles subject to a Gaussian forcing, 
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viscous friction and inelastic collisions. 

The scheme consists of the following ingredients: 



i time is discretized, i.e. t = n ■ dt 

ii update all the velocities to simulate the random forcing and the viscous damping: 



v?(t + dt)=v?(t)e * + yT b (l-e ~)W{t) (29) 

where W(t) is a normally distributed deviate with zero mean and unit variance. 

iii Choose randomly N^- pairs of velocities and update each of them with the collision 
rule (|7[). In this way a mean collision time r c per particle is guaranteed. 

iv Change the time counter n and restart from ii. 

In other words, at every step each particle experiences a Gaussian kick thus receiving 
energy from the bath, whereas it dissipates energy by collision and by damping. For example, 
by choosing dt = 1, mi/Y = t^i = 200 and r c = 25, we obtain that each particle in 
the average experiences 25 Gaussian kicks between two successive collisions and that the 
resulting average kinetic energy is stationary. In order to compare our numerical simulations 
with the theoretical predictions we fixed the temperature of the bath to be TJ, = 1, i.e. 
chosen D = T. The results of such simulations are presented in fig. [| and show a very good 
agreement between the theory and the simulation. 

On the contrary, the agreement between the Sonine expansion and the simulation is not 
completely satisfactory when we consider a system subject to a white noise acceleration, 
but without viscous friction, a driving proposed by some authors |]22j] , [23]. This can be 



considered as the limit r — > 0, T b — > oo, keeping constant D = T b T, in the model defined 
by eqs. (|||): note that in this case the elastic limit r — > 1 cannot be performed without 
taking also the limit D — > as discussed at the beginning, in order to avoid a divergence of 
kinetic energy. For the sake of simplicity, we simulated a one component system (p = 1 and 
m = l) with vanishing viscosity T = 0, but D = 0.0008, r c = 250 and r = 0.5. Such a choice 
yields a granular temperature T g = 16/15, as predicted by our formula (0). Notice that in 
this case the heath bath temperature diverges and the gas does not have a proper elastic 
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limit, since all moments diverge when 7 — > 0. We observed that the tails of the velocity 
distribution function are strongly non Gaussian. These decay as a simple exponential as 
predicted by our simple analysis of the previous section. In fig. (^) we report our simulation 
results against the Sonine approximation. We observe that the theoretical estimate, in spite 
of incorporating the exact values of the first eight moments, deviates from the numerical 
data in the large velocity region. In particular the Sonine expansion can only give Gaussian 
tails, whereas the simulation indicates a slower exponential decay. The reason for such a 
discrepancy is to be ascribed to the slow convergence of the expansion when T = 0. 

.,<> 



o Simulation (t=0) 
□ Simulation (t=3000) 

— Theory (steady value) 
■ ■ ■ ■ Gaussian 

- - Exponential Tail e c 




FIG. 10: Rescaled distribution functions for a one component system (p = 1) with vanishing 
viscosity (r = 0), D = 0.0008 , r = 0.5 and t c = 250. We show the initial Gaussian distribution 
(t = 0) and the asymptotic (t = 3000) stationary distribution. Notice the presence of high velocity 
tails. For the sake of comparison we report the theoretical estimate of the distribution obtained 
by means of the Sonine expansion. 



VI. DISCUSSION AND CONCLUSIONS 

To summarize we have studied the behavior of a model, which perhaps represent the sim- 
plest description of a driven inelastic gas mixture, namely an assembly of two types of scalar 
pseudo-Maxwell molecules subject to a stochastic forcing. We have obtained the velocity 
distributions for arbitrary values of the inelasticity, of the composition and of the masses 
by solving the associated Boltzmann equation by means of a controlled approximation, the 
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moment expansion. The distributions were obtained by computing exactly the moments 
up to the eighth order and then imposing that the corrections to the Maxwell distribution 
stemming from the inelasticity are Gaussians multiplied by a linear combination of Sonine 
polynomials with amplitudes determined self-consistently. 

The model predicts a steady two temperature behavior which is in qualitative agreement 
with existing experimental results. The granular temperatures can be obtained by very 
simple algebraic manipulations for arbitrary values of the control parameters. 

By numerical simulations we demonstrated that the velocity distributions are well de- 
scribed by our series representation in the case of systems in contact with a bath at finite 
temperature T&, whilst the series expansion breaks down in the case of systems in contact 
with bath at infinite temperature, i.e. with zero viscosity. 

What can be learned from such a simple model of granular mixture? Besides obtain- 
ing a global picture of the behavior of the system with a minimal numerical effort both in 
the cooling and in the driven case, the model displays the novel feature of two different 
distribution functions, which remain different even after rescaling by the associated partial 
granular temperatures. Of course the detailed form of the probability velocity distributions 
are strictly model dependent, i.e. depend on the assumption of a constant collision rate in- 
herent in Maxwell models. Finally, the vectorial character of the velocities could be included 
at the cost of a moderate additional effort. A more interesting and difficult problem would 
be that of including in the mixture case a collision frequency proportional to an appropriate 
function of the kinetic granular temperatures, generalizing the work of Cercignani . 

Finally we might ask the general question of the meaning of granular temperature. Our 
findings seem to indicate that it is still the main statistical indicator of the model gran- 
ular system we studied. However, with respect to the temperature of a perfectly elastic 
system it fails to satisfy a very basic requirement which is known as the zeroth principle of 
thermodynamics. 
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APPENDIX A 



In the present Appendix we sketch the derivation of the various moments of the distri- 



bution functions. By equating the equal powers of k in the master equation (11) we obtain 
a set linear of coupled equations for the moments. The method of solution is iterative, be- 
cause the higher moments depend on the lower moments. Thus, for instance, to evaluate the 
lowest order moments of the distribution functions we must solve the following equations 
for the steady state value of the fourth moments: 

(4? - 4r 1 r c )/ii 1) + dgVS 2) + aii(/4 1} ) 2 + ai2/4 1} /4 2) + UD^Jp = (Ala) 
(4? - 4r 2 r c )/4 2 } + + a 22 (/4 2) ) 2 + a 21 ^ ] fif + 12ZVc/4 2) = (Alb) 

In turn, the sixth moments are obtained by solving: 

(4? - 6riT c )/4 1} + 4?/4 2) + hmPfiP + 6i2/i?Vi 2) + b'uvf^P + WDxtjiP = (A2a) 
(4? - 6r 2 r c )/4 2) + 45V? 5 + & 22 4 2) /4 2) + fciA^Vf + MnvPuP + = (A2b) 

Finally the eight moments are the solutions of: 



(d n - »1 it c J/x 8 +a 12 /i 8 + cn/i 2 /i 6 +c u ^ 4 J 

+ r i.W,/ 2 ) 4V // (1) // (2) + r" /; (1) // (2) 4- ^RD r 

l«22 -0 J - 2 r cjA i 8 + "21 ^8 + c 22/^ 2 ^6 +^2^4 J 

+ C2l/^ 2 ^6 + c 2l/ i 4 ^4 + c 2 lA*6 ^2 + OOJJ2T c fl B 

where the general form of the coefficients dij is given by: 



(A3a) 



(A3b) 



4? = -l+pfrS + (1 -7n)1 + (1-P)[7i 2 ] n (A4a) 

4™ ) = (l-p)[(l-7i2)] n (A4b) 

4? = -1 + (1 - P) [7 2 n 2 + (1 - 7 22 )1 + P[72i] n (A4c) 

4i ) =p[(l-7 2 i)] n (A4d) 



and the coefficients are given by: 
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an = 6p[ 7 n(l-7ii)] 2 (A5a) 

a 2 2 = 6(l-p)[ 722 (l- 7 22 )] 2 (A5b) 

a 12 = 6(l-p)[(l-7i 2 )] 2 [7i 2 )] 2 (A5c) 

a 2 i = 6p[(l- 721 )] 2 [ 721 )] 2 (A5d) 



Finally, the coefficients % are given by: 



611 = 
b 2 2 = 

612 = 
bu = 
&21 = 
^1 = 



15p7i 2 i(l -7ii) 2 [7ii + (l-7n) 2 ] 

15(1 - p) 72 2 2 (l - 7 22 ) 2 [7 22 + (1 - 7 22 ) 2 ] 

15(l-p)[(l-7i 2 )] 4 [7i 2 )] 2 

15(l-p)[(l-7i 2 )] 2 [7i 2 )] 4 

15p[(l-7 2 i)] 2 [7 2 i)] 4 

15p[(l-7 2 i)] 4 [7 2 i)] 2 



(A6a) 
(A6b) 
(A6c) 
(A6d) 
(A6e) 
(A6f) 



and Cij are 



c n = 28p7i 2 i(l - 7n) 2 [7 4 i + (1 - 7n) 4 ] (A7a) 

c'n = 70p7n(l-7ii) 4 (A7b) 

c 12 = 28(l-p)[(l-7i 2 )] 6 [7i 2 )] 2 (A7c) 

c' 12 = 70(l-p)[(l-7i 2 )] 4 [7i 2 )] 4 (A7d) 

c'; 2 = 28(l-p)[(l-7i 2 )] 2 [7i 2 )] 6 (A7e) 

c 22 = 28(1 - p) 72 2 2 (l - 7 22 ) 2 [7 2 4 2 + (1 - 7 22 ) 4 ] (A7f) 

c 22 = 70(l-p) 72 4 2 (l- 7 22 ) 4 (A7g) 

c 21 = 28p[(l-7 2 i)] 6 [7 2 i)] 2 (A7h) 

4i = 70p[(l-7 2 i)] 4 [7 2 i)] 4 (A7i) 

c 2 ' 1 = 28p[(l-7 2 i)] 2 [7 2 i)] 6 (A7j) 



It is useful to consider the behavior of the moments in the one component case. The 
major simplicity of the resulting formulae allows us to obtain explicit expressions: 
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Dt c 

H2 = p : r A8a 

rr c + 7(l -7) 

^ 4 = 4 rr c + l-7 4 -(l-7) 4 (A8b) 



30Pr c/M + 15 7 2 (1 - 7 ) 2 (7 2 + (1 - ifmH , AH , 

^ = 6rr c+ l-7 6 -(l-7) 6 (A8C) 



56Dt c ^ 6 + 28 7 2 (1 - 7 ) 2 (7 4 + (1 - 7) 4 We + 70 7 4 (1 - jM rAS ,x 
^ = 8r c r + l-7 8 -(l-7) 8 (AM) 
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